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1.  Introduction 

The  principal  investigator  developed  a  fast  and  computationally 
efficient  implementation  of  a  new  method  for  the  computation  of  steady 
flows  of  viscoelastic  fluids  with  single-integral  constitutive  equations. 
The  P.  I.  and  coworkers  had  developed  a  new  kind  of  finite  element  method 
to  solve  such  problems  {1}  [1,2,3]2.  It  had  been  implemented,  tested,  and 
used  to  solve  problems  of  physical  and  practical  interest.  The  method  was 
computationally  costly,  and  before  it  could  take  its  place  as  a  reliable 
scientific  and  engineering  tool,  the  P.  I.  believed  its  computational  cost 
had  to  be  significantly  reduced,  its  robustness  and  flexibility  improved, 
and  some  measure  of  its  accuracy  in  benchmark  problems  quantified.  The 
aim  of  the  research  was  to  develop  an  optimized  algorithm  that  would  run 
on  the  order  of  a  factor  of  two  faster  on  scalar-architecture  machines 
than  did  the  pilot  algorithm,  and  would  take  full  advantage  of 
vectorization  on  machines  with  that  capability. 

During  the  funded  research,  progress  was  made  in  many  areas,  most 
prominent  were:  improving  the  efficiency  of  the  central  core  of  the 
algorithm,  the  stress  calculator,  and  in  a  understanding  of  the  nature  of 
numerical  difficulties  which  arise  with  some  constitutive  equations  when 
the  stress  calculator  is  used  to  solve  steady  boundary-value  problems.  It 
became  ever  more  apparent  that  the  two  issues  can  be  separated.  It  is 
expected  that  that  the  major  contribution  of  the  completed  research  will 
have  been  the  development  of  a  robust  and  efficient  stress  calculator 
based  on  linear  crossed  triangles.  Such  a  stress  calculator  has  merit  in  its 
own  right  and  can  be  employed  with  a  variety  of  constitutive  equations 
and  with  a  variety  of  numerical  schemes  for  the  solution  of  the  equations 

of  steady  motion  of  the  fluid.  One  such  scheme  is  the  one  the  P.  I.  has  been 

developing,  which  uses  the  crossed-triangle  elements  themselves  to 
develop  a  pseudotime/Galerkin  method.  This  scheme  works  well  with  some 
constitutive  equations:  Deborah  numbers  (non-dimensional  flow  rates)  in 
excess  of  25  can  be  achieved  with  a  Curtiss-Bird  model  {1}  [2-5], 
Deborah  numbers  on  the  order  of  5  with  a  Wagner  [5]/Johnson-Segalman 
{1}  [5,6]  fluid,  on  the  order  of  3  with  an  Oldroyd-B  fluid  [5],  but  only  order 

1  with  a  Maxwell  fluid  [5],  where  the  “High  Weissenberg  Number  Problem” 

seems  most  severe.  A  major  finding  of  the  current  research  is  that  the 
“particle  tracking”  method  of  the  stress  calculator  that  has  been 
developed  can  also  serve  as  the  computational  core  of  new  transient 

'’Numbers  in  braces  refer  to  the  publication  list  of  Section  8.  Those 

in  brackets  refer  to  the  Additional  Reference  Section  10. 
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algorithms  for  non-Newtonian  flows  {10}  with  differential  constitutive 
relations.  Implementation  of  these  techniques  and  their  extension  to 
steady  flows  of  fluids  with  integral  constitutive  equations  is  the  subject 
of  continuing  AFOSR-sponsored  research  (see  also  Section  7.,  below).  It  is 
hoped  that  these  techniques  will  extend  the  range  of  attainable  Deborah 
number  for  those  constitutive  equations  that  currently  have  a  limited 
range. 

2.  Background  and  Problem  Formulation 

The  problem  formulation  and  relevant  equations  are  more  completely 
described  in  Refs.  {1}  [1,2,3];  a  brief  summary  is  given  here.  The 
nondimensionalization  employed  here  is  that  of  Ref.  {11}.  The  complete 
problem  is  typically  to  solve  the  following  set  of  coupled  integro-PDEs, 
subject  to  appropriate  boundary  and  initial  conditions: 

a[-fy-  +  (v-  V)v]  =  V-S 

.  S  =  -  pl  +  2ee  +  t  (2.1) 

T=  /'  (E(E^ - 1) m(f  -  t')dt' 


where  v  is  the  velocity  field,  S  is  the  total  (physical)  stress,  p  the 
hydrostatic  pressure,  e  the  usual  strain-rate  tensor,  and  x  the  non- 
Newtonian  extra  stress.  Ef  is  a  strain  tensor  discussed  in  more  detail 
shortly;  a  is  a  non-dimensional  ratio  of  Reynolds  number  to  Deborah 
number,  and  e  is  the  ratio  of  Newtonian  to  polymer  zero-shear  viscosity 
{5,6,11}.  The  strain-dependence  in  the  history  integral  in  the  last  group  of 
equations  in  Eqs.  (2.1)  is  typical  of  several  kinds  of  integral  constitutive 
equations,  for  example  Johnson-Segalman  {1}  [5,6],  but  the  techniques 
described  here  apply  to  integrands  with  a  much  more  complicated 
functional  dependence  on  E( ,  such  as  is  found,  for  example,  in  the 

Curtiss-Bird  model  {1}  [2-5].  The  function  m(t-t ' )  is  a  function  of  the 
present  time,  t,  and  historical  time,  t ' .  A  typical  memory  function  is  a 
single  exponential, 

m(t  -  t ')  =  e~(,'r)  (2.2) 

where,  following  Refs.  {5,6,9,11},  it  has  been  assumed  that  time  has  been 
scaled  by  the  dominant  relaxation  time.  Though  there  are  many  other 
(perhaps  more  physically  realistic)  possibilities  [5],  Eq.  (2.2)  is  employed 
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in  many  currently  popular  models  and  serves  as  a  representative  example. 

The  key  numerical  difficulty  addressed  in  the  current  research  is  the 
efficient  and  accurate  evaluation  of  the  history  integral.  This  consists  of 
two  separate  problems:  first,  the  integrand  must  be  evaluated,  and  second, 
the  integral  must  be  computed  from  sequence  of  such  evaluations  in 
historical  time,  t'.  As  for  the  first  problem,  the  integrand  is  a  function 
of  the  strain  measured  in  the  fluid  configuration  at  t ' ,  relative  to  the 
configuration  at  t  as  a  reference  state.  This  is  implied  by  an  evolution 
equation  for  the  strain  tensor,  solved  in  reverse  time  along  particle  paths 
emanating  from  the  stress  evaluation  points;  a  common  example  is  given 
by  the  evolution  equation  associated  with  the  ‘non-affine’  deformation 
gradient  of  the  Johnson-Segalman  model  {1}  [5,6], 

)  =  ~  E,(f  ')(ae  -  co)(f ')  (2  3) 

>(0  =  1 


where  co  is  the  spin  tensor,  and  a  is  a  parameter  of  the  model;  when  a=1 
and  Eq.  (2.1)  is  used,  the  Maxwell  model  results  {1}  [5].  To  compute 
particle  pathlines  in  a  time  sequence  of  unsteady  velocity  fields  is 
possible  in  principle,  but  would  require  the  storage  of  a  great  number  of 
velocity  fields.  When  steady  flow  is  assumed  a  priori,  the  problem  is 
greatly  simplified  but  still  formidable.  That  is  the  focus  of  the  current 
research.  In  the  scheme  described  in  Section  4.,  the  time  derivative  in  Eqs. 
(2.3)  is  retained  for  dynamic  relaxation  even  when  steady  flow  is 
assumed. 

The  second  problem  of  evaluating  the  history  integral  is  purely 
numerical  and  will  be  dealt  with  in  Section  4. 

3.  Physical  Problems 

The  techniques  developed  in  the  current  research  apply  directly  to 
isothermal,  planar  flow  problems.  The  P.  I.'s  coworkers  at  Illinois 
Institute  of  Technology  (see  Consultative  and  Advisory  Subsection  of 
Interactions  Section  9.)  have  developed  methodology  for  extending  the 
stress  calculator  techniques  to  axisymmetrical  and  non-isothermal  flows 
[7-9];  implementation  of  these  techniques  in  the  Fast  Algorithm  is  part 
of  the  P.  I.’s  continuing  AFOSR-sponsored  research.  The  flow  problems  that 
have  been  solved  in  the  current  research  are  idealization  of  flows 
encountered  in  polymer  processing  or  property  measurement.  They  serve 
as  test  problems  that  are  complicated  enough  to  fully  test  the  techniques, 
but  simple  enough  so  that  there  is  a  qualitative  understanding  of  what  to 
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expect  of  the  steady  flow  regimes. 

The  first  problem  is  that  of  steady  plane  Poiseuille  flow,  illustrated  in 
Fig.  3.1.  This  is  flow  between  two  flat,  infinite  plates,  driven  by  a 
constant  pressure  gradient.  The  primary  interest  in  this  flow  here  is  as 
boundary  data  for  the  truly  two-dimensional  flows.  The  whole  question  of 
well-posedness  of  planar  flow  problems  is  a  delicate  one  and  is  only 
understood  for  steady  flows  with  certain  constitutive  equations  [10].  The 
P.  I.  has  performed  extensive  numerical  experiments  and  found  that 
imposing  mass  conserving  velocity  profiles  at  inlets  and  outlets  that  are 
essentially  “at  infinity"  works  well  in  most  cases.  Further  study  of  these 
matters  is  part  of  the  P.I.’s  NSF-sponsored  research  (see  Support 
Statement  Section  12.).  That  research  also  involves  the  study  of  the 
dynamics  of  shear  flows  leading  to  steady  profiles  such  as  that  of  Fig.  3.1 
{5,6,9,11,13}.  An  important  insight  gained  in  that  research  is  that  the 
steady  flow  calculations  described  here  cannot  be  expected  to  converge  in 
certain  cases  because  of  the  development  of  singularities  in  the  shear 
rate.  This  matter  is  discussed  in  more  detail  below. 
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Figure  3.1.  A  plane  Poiseuille  flow  profile.  This  flow  is  used  as  boundary  data 
in  some  of  the  planar  flows.  Symmetry  is  assumed  in  this  figure,  but  the  full 
profile  is  also  useful  as  boundary  data  in  cases  where  enforced  symmetry  is  not 
appropriate. 
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Fig.  3.2  illustrates  planar  contraction  flow.  The  pictured  domain  is  a 
cross-section  of  a  wide  channel  with  an  abrupt  change  in  height;  only  the 
upper  half  of  the  channel  is  modelled  because  symmetry  of  the  geometry 


inflow 


Figure  3.2.  Abrupt  contraction  flow  test  problem,  showing  a  very  crude  mesh. 

Inflows  and  outflows  are  assumed  to  be  far  enough  from  the  region  of  disturbed 
flow  so  that  plane  Poiseuille  flow  can  be  imposed  (see  Fig.  3.1).  The  inset  shows 
a  square  special  case  of  the  quadrilateral  crossed  triangle  macroelement  used  in 
this  research;  it  is  composed  of  four  constant  strain-rate  triangles  formed  by 
crossing  quadrilateral  diagonals. 

and  resulting  flow  is  assumed.  A  major  reason  for  studying  the  abrupt 
contraction  flow  is  that  the  singularity  in  shear  rate  induced  by  the 
corner  appears  to  be  significant,  though  its  order  is  not  known  for  most 
non-Newtonian  fluids  [1]. 

The  journal  bearing  problem  illustrated  in  Fig.  3.3  is  an  idealization  of 
an  important  lubrication  flow.  In  most  actual  lubrication  applications 
(main  bearings  in  internal  combustion  engines,  for  example),  the 
lubricating  film  is  so  thin  that  the  gap  between  bearing  and  shaft  is  too 
narrow  to  be  modelled  by  the  techniques  described  here;  however,  the 
eccentric  geometry  is  of  interest,  because  for  larger  gaps  it  provides  a 
problem  for  which  spectral  techniques  can  provide  accurate  numerical 
solutions  at  very  high  Deborah  number  [11].  These  solutions  can  be  used  as 
benchmarks  for  fully  discrete  methods;  this  has  led  others  to  develop  new 
numerical  formulations  for  differential  constitutive  equations  that  are 
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very  successful  at  high  Deborah  number  [12].  It  is  hoped  that  similar 
benefits  to  integral  formulations  will  accrue  in  the  future.  One  of  the 
most  important  points  about  the  journal  bearing  problem  is  that  there  are 
no  corners  in  the  domain  and  no  inflows  or  outflows.  The  concentric  case 
‘true  Couette  flow’  is  of  interest  because  the  exact  solution  is  known. 


Figure  3.3.  Domain  for  journal  bearing  problem.  This  highly  idealized 
approximation  to  a  lubrication  flow  is  a  standard  test  problem  in  which  a  wide 
variety  of  gaps  an  eccentricities  can  be  specified.  The  most  numerical  difficulty 
occurs  for  high  eccentricity  and  narrow  gaps. 
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Because  the  eccentric  cylinder  geometry  is  of  such  mathematical  and 
numerical  interest,  experimenters  have  been  led  to  make  flow 
visualization  emparati  incorporating  this  geometry  [13],  and  the  P.  I.  has 
extensive  L  _*ser  Doppler  Velocimetry  data  for  such  flows  at  his  disposal. 
Making  a  comparison  between  computed  and  observed  solutions  is  an 
important  component  of  the  P.  I.’s  continuing  AFOSR-sponsored  research. 


Figure  3.4.  Mesh  for  journal  bearing  problem  in  special  case  when  shaft  and 
bearing  are  concentric;  exact  solution  is  known  for  this  test  problem.  For  the 
crossed  triangle  macroelement  of  Fig.  3.2,  using  an  odd  number  of  elements  in 
each  circumferential  strip  is  essential  to  avoid  ‘checkerboard’  problems. 
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The  problem  to  which  most  effort  in  the  current  research  has  been 
devoted  is  that  of  a  plane  flow  over  a  transverse  slot,  pictured  in  Fig. 
3.5.  If  the  pictured  flow  is  driven  by  motion  of  the  upper  plate,  we  say 
the  flow  nas  “Couette  base  flow”;  if  it  is  driven  by  an  applied  pressure 
gradient  while  both  walls  stay  fixed,  we  say  that  it  has  “(plane) 
Poiseuille  base  flow.”  This  flow  (with  Poiseuille  base  flow)  is  the 
basis  of  a  unique  in-line  measurement  device  for  measuring  viscosity 
anc'  first  normal  stress  differences  of  polymer  solutions  and  melts,  the 
Lodge  Stressmeter™  {1,2,8,12}  [2,3,5,14].  The  device  uses  the  ‘elastic 
hole  pressure’  [2,3,14]  to  predict  the  first  normal  stress  difference. 
The  elastic  hole  pressure  is  the  difference  in  thrust  between  the  top 
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Figure  3.5.  Domain  for  slot  flow  problem;  h,  b,  and  d  are  physical  dimensions 
that  correspond  to  actual  device  parameters.  L  is  assumed  to  be  large  enough  so 
that  plane  Poiseuille  flow  can  be  imposed  at  the  in-and  outflows  to  a  reasonable 
degree  of  accuracy  (see  Fig.  3.1). 

and  bottom  of  the  slot,  minus  the  inertial  contribution  to  the  hole 
pressure  [2,3,14].  The  shear  viscosity  is  measured  by  a  rather  straight¬ 
forward  estimation  of  the  pressure  gradient  [14]  [In  fact,  actual 
devices  have  two  widely  separated  slots  for  this  purpose.],  but  the 
normal  stress  measurement  is  based  on  a  rather  subtle  measurement 
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relation  [15].  A  major  reason  for  studying  this  flow  is  the 
“Higashitani-Pritchard-Baird-Lodge  (HPBL)  Paradox”:  The  derivation  of 
the  HPBL  measurement  relation  in  Ref.  [15]  can  only  be  approximately 
correct  in  Couette  base  flow  {17};  in  the  Poiseuille  base  flow,  the 
derivation  leaves  out  terms  that  are  0(1)  with  respect  to  the  measured 
quantities  {12}  [16],  yet  the  measurement  relation  appears  to  be 
correct.  This  paradox  was  uncovered  by  the  P.  I.  and  his  Research 
Assistant  during  the  funding  period  (and  independently  by  R.  Srinivasan 
[16]);  it  now  appears  to  be  resolved  (as  outlined  in  the  Status  of 
Research  Section  5.).  This  aspect  of  the  problem  is  a  major  component 
of  the  thesis  work  described  in  the  Publications  Section  8  (Ref.  {17}). 


Figure  3.6.  Typical  mesh  for  slot  flow  problem.  Grid  is  finest  in  the  slot 
mouth  area,  where  experience  shows  that  the  curvature  of  the  streamlines  in 
the  channel  and  at  the  top  of  the  vortex  in  the  slot  are  crucial  to  accurate 
determination  of  the  hole  pressure. 

During  the  funding  period  other  aspects  of  slot  flow  have  been  studied, 
particularly  the  excess  pressure  rise  in  the  base  Couette  or  Poiseuille 
flow  due  to  the  slot  disturbance.  The  slot  flow  problem  is  also  interesting 
for  the  qualitative  features  observed.  One  such  feature  is  the  appearance 
of  vortices  in  the  slot  itself;  these  have  been  the  subject  of  extensive 
flow  visualization  study  [17], 
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4.  Tools  and  Techniques 

Equations  (2.1)  are  solved  by  a  standard  Galerkin  formulation,  based  on 
the  crossed  triangle  macroelement  illustrated  in  Fig.  3.2  [1,2,3].  The 
history  integral  is  evaluated  by  a  weighted  quadrature  formula  that 
includes  the  memory  function  in  the  weights.  When  the  memory  function  is 
the  single  exponential  of  Eq.  (2.2),  the  quadrature  formula  is  a  Laguerre 
formula,  otherwise  it  is  some  appropriate  generalization  of  a  Laguerre 
formula  [2].  To  evaluate  the  integrand  of  the  history  integral  of  Eqs.  (2.1) 


2 


Figure  4.1.  A  constant  strain-rate  triangle  showing  a  conic  section 
streamline/particle  path.  Arrows  indicate  flow  direction.  Filled  circle  is  the 
location  of  the  centroid/stress  evaluation  point  in  the  first  element  and  the 
location  at  which  the  generalized  Laguerre  point  is  found  in  succeeding  elements 

(not  necessarily  the  centroid);  t  c  is  0  for  the  stress  evaluation  element  and  is 

a  generalized  Laguerre  time  otherwise;  t  'e  and  t  '0  are  element  boundary 
crossing  times. 

“particle  tracking”  is  used.  This  technique  is  described  in  detail  in  [1,2,3], 
and  is  illustrated  in  Fig.  4.1.  Tracking  is  performed  from  centroid  of 

triangle  where  stress  is  evaluated  with  historical  time  t'c  =  0,  following 
path  A  backwards  along  the  upstream  portion  of  the  streamline,  emanating 
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from  the  centroid  at  the  filled  circle.  In  succeeding  elements  the  paths 
are  like  the  union  of  portions  A  and  B,  with  boundary  crossing  times  t ' 
and  t  '.  This  tracking  is  continued  back  into  the  past  until  such  time  as  a 
temporal  generalized  Laguerre  point  is  reached  at  time  t '  at  the  spatial 

location  of  the  filled  circle  (now  not  necessarily  the  element  centroid). 
The  integrand  is  evaluated,  and  the  process  continues  until  the  last 
integration  point  is  encountered. 

The  crossed  triangle  macroelement  is  ideally  suited  to  this  procedure 
because  it  has  excellent  mass  conservation  properties  and  constant 
velocity  gradients  on  each  triangle  [2,3].  In  this  case,  one  finds  that  Eqs. 
(2.3)  have  an  exact  solution  on  each  triangle,  at  any  historical  time  t ' 
between  t  '  and  t  ' .  This  solution  is  given  by 


E,(f  0  -  (cosh[ SO' -t  ;)])E,( f ;)  -  Sinh[^'-';)]E ,( f '.)  A 


(4.1) 


where  A  is  the  coefficient  matrix  in  Eqs.  (2.3)  and  5  is  its  eigenvalue 
(with  a  plus  sign).  Using  the  initial  condition  of  Eqs.  (2.3)  and 
accumulating  the  strain  multiplicatively  according  to  Eq.  (4.1),  the 
integrand  can  be  calculated  analytically  along  the  streamline.  Two 
essential  features  of  these  calculations  are  [2,3]:  1)  The  streamlines  are 
conic  sections  on  each  element  and  C-j  at  element  boundaries,  and  2)  there 

is  an  analytic  relation  between  any  two  points  x-j  and  X2  on  a  streamline 
and  the  transit  time  between  these  points,  given  by  the  drift  function, 


t'-  t'2=  w(Xl)-  w(x2) 


(4.2) 


The  calculations  just  described  allow  evaluation  of  the  history  integral, 
thus  the  total  stress  at  any  triangle  centroid  in  the  domain.  This 
procedure  is  embedded  in  an  implicit-explicit  pseudotime/Galerkin 
scheme,  generalizing  that  of  Hughes,  Liu,  and  Brooks  [18].  In  the  following 
equations,  the  time  level  is  denoted  by  subscripts  n  or  n+1  and  inner 
iterations  by  superscripts  (0),  (/),  or  (41).  Absence  of  superscripts 
indicates  fully  converged  inner  iteration. 
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Ma^,+Cv,.,+  P<V„J  +  Q(V„J = F„„ 

P(v)  =  aj  Nr(v  ■  Vv)dV 

n 

Q(v)  =  J  BTSdV  -  Cv 

a 

vl°i.  =  v»  +  (1- 1')4,a»  (4.3) 

v1::: =v(::,-j'’{[m+ 

+  rA(P(v^,)  -Mv™  +  yAfF^,} 

an.i  =  (v»..-v'l)/5'Af 
J=M+yAtC  +  (yA(-§£)  +  (yAf-|S) 

opt  opt 


M  and  C  are  the  usual  finite  element  mass  and  dissipation/penalty  {1,2,3} 
matrices,  corresponding  to  the  nondimensionalization  of  Eqs.  (2.1)  or  in 
analogous  dimensional  form.  P  is  the  nonlinear  contribution  to  the 
material  derivative,  and  Q(v)  is  the  non-Newtonian  contribution  to  the 
total  stress,  with  S  computed  from  the  indicated  velocity  field  by 
particle  tracking.  N  is  the  matrix  of  shape  functions  and  B  the  usual  finite 
element  ‘B-matrix’  {18}.  J  is  an  approximation  to  the  Jacobian  of  the 
residual  expression  in  the  braces  following  its  inverse,  and  y  is  parameter 
of  the  integration  scheme  usually  chosen  to  be  between  0.5  and  1.  The  ‘opt’ 
subscripts  on  the  second  two  terms  in  the  definition  of  J  indicate  that 
this  Jacobian  information  is  optional.  In  the  current  algorithm  these 

terms  are  omitted,  in  order  to  preserve  symmetry  and  simplicity  of  the 
matrices  involved.  Two  compensatory  techniques  are  employed:  1)  An 
inverse  Broyden  quasi-Newton  scheme  is  used  that  has  the  effect  of 
approximating  these  terms  as  the  inner  iterations  proceed,  and  2)  an 
upwinding  scheme  generalizing  that  of  Ref.  [18]  to  constant  strain-rate 
triangles  is  used  on  the  inertial  terms.  Though  the  upwinding  scheme  is 
not  particularly  sophisticated,  it  seems  to  work  very  well  at  the  lower 
Reynolds  numbers  characteristic  of  the  very  viscous  flows  studied  in  this 
research. 

When  Q*0,  the  integration  scheme  of  Eqs.  (4.3)  is  a  ‘pseudotime 
scheme’  because,  though  the  time  derivative  in  the  inertial  term  is 

retained  and  differenced,  the  time-dependence  of  Q  is  not  properly  taken 
into  account.  The  stress  is  calculated  assuming  that  the  current  velocity 
field  takes  the  place  of  all  velocity  fields  in  the  past,  so  that  these  need 

not  be  stored.  This  assumption  is  only  valid  when  steady  flow  is 
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to  the  steady  state.  This  feature  must  be  viewed  as  a  limitation  of  the 
approach,  but  the  alternative  of  fully  dynamic  simulation  with  integral 
models  seems  out  of  reach  at  present.  The  hope  is  that  the  algorithmic 
damping  of  the  pseudotime  scheme  will  lead  to  a  more  rapid  convergence 
at  higher  Deborah  numbers  than  can  be  achieved  with  direct  steady 
iteration.  It  should  be  noted  that  direct  steady  iteration  can  be  obtained 
from  Eqs.  (4.3)  as  Af->°°  by  taking  y=1. 

The  key  issue  in  the  algorithm  of  Eqs.  (4.3)  is  the  efficient  evaluation 
of  the  extra  stress,  for  each  of  many  inner  iterates.  This  is  a 
computationally  burdensome  calculation,  since  many  generalized  Laguerre 
points  are  required  to  obtain  accurate  stresses;  this  implies  that 
thousands  of  element-level  tracking  calculations  must  be  carried  out  for 
each  iteration. 


5.  Status  of  the  Research 

During  the  reporting  period,  research  progress  has  been  made  in  several 
areas: 

1.  Vectorization  of  the  fast  algorithm:  The  vectorization  of  the  P.  I.’s  code 
for  two  supercomputers  has  been  completed,  and  the  P.  I.  and  his 
Research  Assistant  are  in  the  process  of  preparing  a  version  of  the 
code  for  the  SDSC  Cray  XMP/48.  The  resulting  improvements  are 
detailed  in  Section  6.  As  expected,  since  particle  tracking  involves  so 
many  element  level  calculations,  vectorization  did  not  make  as  much  a 
relative  difference  for  non-Newtonian  flow  as  it  did  for  Newtonian 
(Navier-Stokes)  flow.  Other  enhancements  described  in  the  next  sub¬ 
section  proved  to  be  more  effective  in  the  non-Newtonian  case. 

2.  Code  enhancement  and  streamlining:  In  current  supercomputing 
environments  core  storage  capacity  is  enormous  with  respect  to  the 
needs  of  the  P.  I.'s  code;  this  is  also  true  in  virtual  memory 
environments.  It  has  proved  extremely  cost-effective  to  avoid 
regeneration  of  element-level  quantities  wherever  possible,  in  favor  of 
arrays  that  are  permanent  or  redefined  only  at  the  outset  of  each 
iteration.  Unlike  standard  finite  element  procedures,  those  for 
memory  fluids  must  repeatedly  use  information  which  is  nonlinearly 
dependent  on  the  current  the  velocity  field  at  the  current  time  step.  For 
standard  finite  element  procedures,  the  choice  is  to  generate  the 
information  as  needed,  discarding  the  result  after  it  is  used.  Here,  it 
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pays  to  generate  element  data  such  as  strain  rate  and  save  it  in  a  large 
array  at  the  outset  of  each  inner  iteration  at  each  pseudotime  step. 
Since  each  element  can  contain  the  pathlines  of  many  material  points 
required  to  compute  the  new  stress  field,  numerous  regenerations  of 
the  data  from  nodal  values  can  be  avoided.  For  the  P.  I.  's  code  the 
saving  attributable  to  this  trade-off  was  found  to  be  20%  to  compute  a 
solution  in  abrupt  contraction  flow  for  a  Johnson-Segalman/Wagner 
fluid  at  Deborah  number  5. 


3.  Adaptive  memory  quadrature:  This  is  the  technique  of  determining  the 
number  of  generalized  Laguerre  points  for  each  triangle,  based  on  the 
maximum  shear-rate  encountered  on  the  particle  path  emanating  from 
the  element  in  the  previous  iteration.  A  simple  thresholding  scheme 
was  used  in  which  all  triangles  with  relatively  high  shear  rates  use  the 
maximum  available  number  of  quadrature  points,  but  triangles  with 
significantly  lower  shear  rates  use  a  reduced  number.  It  has  been  found 
that  a  substantial  reduction  in  running  time  can  be  obtained  with  the 
Maxwell  fluid  model  from  the  adaptive  quadrature  alone.  The  saving  is 
less  for  more  complicated  constitutive  equations  but  still  is 
significant.  For  the  Johnson-Segalman/Wagner  fluid  in  abrupt 
contraction  flow  at  Deborah  number  5,  the  saving  due  to  adaptive 
quadrature  was  found  to  be  about  20%.  Note  that  the  combined  savings 
of  items  2.  and  3.  amount  to  a  reduction  of  40%  in  CPU  cost  in  the 
Johnson-Segalman/Wagner  case  and  much  more  in  the  Maxwell  fluid 
case. 


4.  Model  problems  for  the  High  Weissenberq  Number  Problem  (HWNP);  This 
term  refers  to  the  pathologies  which  cause  difficulty  in  converging  to 
solutions  with  significant  non-Newtonian  effects  when  many 
constitutive  equations  are  employed  [1].  The  P.  I.  has  succeeded  in 
developing  a  simple,  analytically  tractable  model  which  shows  possible 
sources  of  the  HWNP  {13}.  In  the  model  problem,  the  combined  equations 
of  steady  motion  and  stress  lose  their  ellipticity  with  increasing  flow 
rate.  This  can  be  treated  in  1-D  by  techniques  that  are  closely  related 
to  optimal  upwinding  techniques  for  the  Navier-Stokes  equations.  In 
two  space  dimensions  the  problem  is  more  difficult:  Added  viscosity 
must  be  anisotropic;  previous  attempts  to  add  ‘Newtonian’  viscosity 
seem  to  have  failed  for  this  reason.  It  is  also  unclear  how  to  define  a 
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‘Peclet-like’  number  for  planar  flows;  anisotropy  indicates  that  there 
should  be  at  least  two  of  them. 

Further  attempts  to  extend  the  concept  of  anisotropic  viscosity 
were  deferred  because  of  two  interesting  discoveries:  First,  in  the 
analysis  of  the  equations  of  motion  and  stress  for  a  Johnson-Segalman 
model  with  a  very  small  amount  of  ‘Newtonian  viscosity,’  it  was  found 
that  the  full  system  had  an  unsteady  transition  to  asymptotically 
steady  ‘spurt’  solutions  with  ‘singular  surfaces’  {5,6,9,11,14}  [19]. 
These  solutions  are  believed  to  be  physically  important  and  can  be 
related  to  the  phenomenon  of  ‘melt  fracture’  [20].  The  added  viscosity 
techniques  would  artificially  suppress  such  phenomena,  because  their 
signature  is  precisely  the  loss  of  ellipticity  in  the  steady  equations. 
Second,  for  spurt  solutions  in  shear-flow,  numerical  methods 
suggested  themselves  which  require  no  anisotropic  viscosity,  which 
are  fully  dynamic,  and  are-  robust  to  the  development  of  singular 
surfaces.  The  extension  of  these  methods  to  2-D  flows  is  the  central 
focus  of  this  continuing  AFOSR-sponsored  research.  The  key  insight  is 
that  the  loss  of  ellipticity  of  the  steady  equations  is  only  fatal  to 
methods  which  ‘parabolize’  the  steady  flow  problem,  either  by  direct 
stationary  iteration  —  which  is  analogous  to  a  pseudotime  scheme  for 
parabolic  equations  —  or  by  the  overt  introduction  of  such  a 
pseudotime.  Ellipticity  loss  is  not  fatal  to  schemes  which  respect 
hyperbolic  character,  and  the  P.  I.  proposes  that  the  new,  efficient 
stress  calculator  can  be  imbedded  in  such  schemes. 


5.  Exactly  incompressible  projection:  This  is  a  technique  for  the  removal 
of  residual  compressibility  from  finite  element  interpolations  to 
velocity  fields.  It  was  developed  by  E.  T.  Olsen  (Dept.  Mathematics, 
l.l.T.)  as  part  of  the  P.  I.'s  joint  research  with  the  l.l.T.  group  (see 
‘Interactions’  Section  9.  following).  The  idea  is  to  project  the  slightly 
compressible  field  onto  the  weakly  incompressible  fields  (which,  for 
the  crossed  triangle  element,  are  exactly  incompressible)  in  some 
convenient  norm.  The  normal  equations  are  solved  by  conjugate- 
gradients  without  assembly  of  the  system  matrix.  Since  the  correction 
is  typically  very  small,  the  process  is  rapidly  convergent.  The  need  for 
projection  arises  from  the  fact  that  the  ‘tracking’  procedure  assumes 
exact  incompressibility.  When  the  slight  compressibility  arises  from, 
the  penalty  used  to  give  pressure  approximations,  it  is  negligible  for 
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the  purpose  of  tracking,  except  in  some  extreme  and  rarely  occurring 
cases.  But  when  the  stress  calculator  is  used  on  experimentally 
observed  velocity  fields,  experimental  and  interpolation  error  could 
lead  to  significant  artificial  compressibility  which  must  be  eliminated 
before  the  stress  calculator  can  function  reliably. 

6.  Development  of  multimesh  extrapolation  procedures:  In  the  early  phases 
of  the  current  research,  the  P.  I.  was  puzzled  by  the  fact  that  numerical 
solutions  failed  to  confirm  the  HPBL  measurement  relation  (see  Section 
3.)  [2,3,5].  M.  F.  Webster  (formerly  Reading,  U.  K.,  now  University 
College,  Swansea)  reported  good  agreement  between  finite  difference 
solutions  and  the  HPBL  relation.  At  the  P.  I.’s  encouraging,  Webster 
performed  grid  refinement  studies  and  found  that  the  results  changed 
and  there  was  poorer  agreement  with  the  HPBL  relation  on  finer  grids. 
The  P.  I.  and  Webster  collaborated  in  the  development  of  Richardson 
extrapolation  techniques  to  investigate  the  behavior  of  our  respective 
numerical  schemes  {2}.  We  found  that  for  the  Maxwell  model,  for  which 
both  of  our  numerical  methods  could  be  used,  exhibits  slow  first  order 
convergence  of  the  modelled  slot  pressure,  using  the  P.  I.’s  finite 
element  method,  and  non-monotonic  second-order  convergence  for  the 
finite  difference  scheme.  Both  methods  appeared  to  be  converging  to 
the  same  hole  pressures  and  related  quantities,  but  would  require  very 
fine  grids  to  get  produce  hole  pressures  in  line  with  the  HPBL 
predictions;  however,  the  extrapolated  results  agreed  well  between  the 
two  methods  and  gave  good  agreement  with  the  HPBL  relation. 
Agreement  with  the  HPBL  relation  for  the  Maxwell  fluid  in  planar  slot 
flow  was  later  confirmed  by  Sugeng,  Phan-Thien,  and  Tanner  [21],  using 
a  boundary  element  method.  With  some  other  constitutive  equations, 
such  as  Curtiss-Bird  [4],  in  which  disagreement  with  the  HPBL  relation 
is  observed,  the  error  indicators  developed  in  conjunction  with  the 
extrapolation  techniques  yield  ambiguous  results,  and  their 
satisfaction  of  the  HPBL  relation  is  still  an  open  question. 


7.  The  solution  of  problems  of  physical  interest:  In  spite  of  the  Deborah 
number  limitations  of  the  current  pseudotime  approach,  significant 
progress  was  made  in  the  understanding  of  several  interesting 
problems.  In  flows  over  transverse  slots  it  was  discovered  that,  while 
the  HPBL  measurement  relation  is  apparently  misapplied  to  pressure- 
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gradient  driven  flows,  its  predictions  are  evidently  correct.  The  thesis 
work  of  M.  Yao  {17}  shows  promise  of  explaining  this  paradox.  We  find 
that  the  derivation  of  the  HPBL  relation  contains  two  steps  that  are 
only  approximately  correct;  with  a  Couette  base  fiow,  for  a  Maxwell  or 
Oldroyd-B  [5]  fluid,  the  approximations  are  accurate,  but  with  plane 
Poiseuille  base  flow  they  are  not;  however  the  error  from  the  two 
questionable  steps  cancels,  leading  to  a  correct  result.  Mr.  Yao  is 
continuing  to  investigate  the  possibility  of  error  cancellation  with 
other  constitutive  equations.  One  of  the  approximate  steps  in  the  HPBL 
analysis  is  the  dropping  of  the  streamwise  gradient  of  the  normal 
stress.  One  might  suppose  that  this  is  a  reasonable  assumption  with 
Couette  base  flow,  but  we  have  found  that  there  can  be  a  significant 
pressure  rise  in  such  flows,  induced  by  the  disturbance  of  the  slot. 
Fortunately,  in  the  derivation  of  the  HPBL  relation,  the  integral 
contribution  from  this  pressure  rise  is  small  with  Couette  base  flow. 
We  have  studied  these  slot-induced  pressure  rises  with  both  Couette 
and  Poiseuille  base  flows,  and  details  of  the  investigation  will  be 
found  Ref.  {17}.  We  have  also  computed  flows  in  abrupt  contraction 

Slot  Flow  Contraction  Flow  Journal  Bearing 

Maxwell  Stress  C.  Maxwell  Stress  C.  Maxwell  Stress  C. 


De 

0.9167 

3.3040 

1.3290 

4.7520 

4.4165 

8.8500 

Iters. 

7 

12 

9 

8 

4 

10 

CPU 

230 

394 

174 

207 

123 

700 

Table  5.1  Summary  of  maximum  Deborah  number  (De)  and  CPU  time  in 
minutes  achieved  by  the  pseudotime/Galerkin  method  in  various  flows.  'Stress 
C.'  refers  to  a  single  integral  model  with  similar  characteristics  to  a  modified 
Phan-Thien/Tanner  model  [5].  Slightly  higher  De  can  probably  be  achieved  in 
slot  flow.  Deborah  numbers  in  excess  of  25  can  be  achieved  with  this  method 
using  the  Curtiss-Bird  [4]  constitutive  equation,  but  the  new,  fast  stress 
calculator  has  not  been  implemented  to  use  that  constitutive  equation  yet.  The 
mesh  for  slot  flow  has  1008  macroelements  {2},  for  contraction  fiow  940,  and 
for  journal  bearing  flow  384  macros.  CPU  times  are  for  a  VAX  11/780  using 
double  precision  arithmetic. 

and  in  planar  idealizations  of  journal  bearings.  In  the  latter  case,  we  are 
involved  in  a  detailed  comparison  of  our  results  with  experimental 
velocity  measurements  (the  thesis  work  of  S.  Burdette  [13]).  Table  1  gives 
a  summary  of  the  maximum  Deborah  numbers  achieved  and  the 
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computational  cost  on  a  Vax  11/780.  A  more  complete  description  of 
these  results  and  others  will  be  found  in  Ref.  {17}. 

6.  Performance  Improvement 

The  central  focus  of  the  funded  research  was  to  improve  the 
performance  of  the  original  pilot  code  and  enhance  its  usefulness  as  a 
research  tool.  The  following  summary  is  excerpted  from  Ref.  {17}.  The 
pilot  code  referred  to  as  FLUCODE  Version  0.0  had  already  been  modified 
and  enhanced  by  the  P.  I.  under  the  current  project  (for  example,  a 
preliminary  version  of  adaptive  memory  quadrature  had  been  implemented) 
and  thus  represents  a  more  advanced  version  of  the  code  than  the  original 
pilot  version.  The  performance  improvement  documented  below  in  addition 
to  the  enhancements  already  made  in  Version  0.0  add  up  to  an  average 
performance  improvement  in  excess  of  the  project  target  of  a  factor  of 
two  improvement  in  running  time.  For  example,  in  the  preparation  of  Ref. 
{2},  a  run  comparable  to  Column  1  of  Table  5.1  was  made  with  the  basic 
pilot  code.  It  required  16  iterations  and  677  CPU  minutes.  According  to 
Table  5.1,  the  enhanced  algorithm  has  reduced  the  number  of  required 
iterations  by  more  than  a  factor  of  two  and  reduced  the  CPU  time  per 
iteration  by  22%  over  the  basic  pilot  version. 

Mr.  Yao  has  also  developed  an  interface  with  FIDAP  [22]  for  pre-and 
postprocessing.  This  fluid  dynamics  analysis  package  is  widely  used  and 
makes  FLUCODE  more  easily  accessible  to  other  researchers.  It  is 
interesting  to  note  that  FLUCODE  and  FIDAP  can  both  solve  Navier-Stokes 
problems,  and  that  the  Navier-Stokes  solver  in  FLUCODE  performs  very 
well  in  comparison  to  FIDAP.  It  should  be  pointed  out,  however,  that  FIDAP 
has  many  general  purpose  features  that  FLUCODE  does  not. 

The  versions  of  FLUCODE  are  designated  as  follows: 

•Version  0.0  -  pilot  code  developed  by  D.S.  Malkus,  et  al. 

•Version  0.1  -  vectorized  code  of  version  0.0; 

•Version  1.0  -  the  earlier  stage  of  the  new  version  (modified 
by  M.  Yao)  without  storing  the  geometrical,  material  and 
iterational  quantities; 

•Version  1.1  -  the  newest  VAX  version  of  FLUCODE  (modified 
by  M.  Yao); 

•Version  1.2  the  newest  CRAY-2  version  of  FLUCODE 
(modified  by  M.  Yao). 

The  test  problem  referred  in  this  section  is  the  4:1  Abrupt  Contraction 
Flow  Simulation  of  Fig.  3.2,  described  in  Section  3.  The  symmetry  of  flow 
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domain  is  assumed  and  three  cases  are  considered: 


•Stokes  solution,  Re=0,  linear  case; 

•Navier-Stokes  solution,  Re=80,  56  explicit  time  steps, 
At=0.02; 

•Modified  Johnson-Segaiman  fluid,  De=0.3, 

‘pseudotime  step  with  predictor-corrector  plus  1  inner 
iteration,  At=1.0, 

‘explicit  pseudotime  steps,  At=1.0. 

Significant  performance  improvement  has  been  achieved  by  the  newest 
version  of  FLUCODE.  Table  6.1  shows  the  performance  improvement  on  the 
VAX  11/780  and  performance  comparison  with  FIDAP,  for  the  Stokes 
solution  case.  The  improvement  of  Version  1.1  is  about  21%  over  Version 
0.0.  The  performance  improvement  for  the  Navier-Stokes  solution  case  can 
be  seen  in  Table  6.2  &  6.3.  There  is  about  25%  CPU  saving  for  the  Navier- 
Stokes  solution  for  Version  1.1  on  the  VAX.  Tables  6.4  to  Table  6.6  show 
the  modified  Johnson-Segaiman  fluid  case.  From  Table  3.5  we  can  see  that 
the  performance  improvement  of  version  1.2  on  the  CRAY-2  is  about  34% 
over  Version  0.1.  This  is,  in  some  part,  due  to  storing  the  geometrical  and 
material  information  at  the  beginning  and  the  intermediate  iteration 
information  at  each  time  step  instead  of  re-calculating  them.  As  a  price, 
the  in-core  storage  needed  in  the  newest  version  is  much  larger  than 
Version  0.0. 


rmi 

FLUCODE  0.0 
FLUCODE  1.0 
FLUCODE  1.1 
FIDfiP  2.0 


CPU  (seconds) 

135.8 

110.5 

107.0 

226.3 


Table  6.1  Performance  Improvement  on  VAX-11/780 
4:1  Abrupt  Contraction  Flow  (940  elements,  1024 
nodes)  Stokes  Solution  (Rc=0). 
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CODE 

Iterations 

CPUfsec.1 

FLUC0DE  0.0 

56 

1646.3 

FLUC0DE  1.1 

56 

1239.1 

FLUC0DE  1.1 

1 

197.0 

FIDflP  2.0 

1 

373.0 

Table  6.2  Performance  Improvement  on  VAX  11/780 
4:1  Abrupt  Contraction  Flow  (940  elements,  1024 
nodes)  Navier-Stokes  Solution  (Re=80). 


CRflY-2  at  U.  of  Minnesota 


ELUXODL 

Comoiler 

Uector 

Iters. 

CPUfsec.) 

0.0 

cft77  0.23 

disabled 

56 

81 

0.1 

cft77  0.23 

disabled 

56 

79 

0 . 1 

cft77  0.23 

enabled 

56 

57 

1.2 

cft77  1.3a 

enabled 

56 

— 

Table  6.3  Performance  Improvement  on  CRAY-2  4:1  Abrupt 
Contraction  Flow  (940  elements,  1024  nodes  )  Navier-Stokes 
Solution  (Re=80). 


UflH-1 1/780  at  CMS  —  U.  of  Wisconsin 


E1U.CQDL 

Compiler 

Iterations 

C£U(sec..) 

0.0 

Fortran 

l 

2881 

1.0 

Fortran 

l 

2407 

1 . 1 

Fortran 

l 

— 

0.0 

Fortran 

3 

4862 

1.0 

Fortran 

3 

3906 

1 . 1 

Fortran 

3 

3348 

Table  6.4  Performance  Improvement  on  VAX-11/780  4:1 
Abrupt  Contraction  Flow  (940  elements,  1024  nodes)  Modified 
Johnson-Segalman  fluids  with  De=0.3. 
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CRflY-2  at  U.  of  Minnesota 


FLUCODE 

Comoiler 

Uector 

Iters. 

CPUfsec.l 

o 

o 

eft 77  0.23 

disabled 

3 

158 

0.1 

eft 77  0.23 

disabled 

3 

157 

0.1 

cft77  0.23 

enabled 

3 

139 

1.2 

cft77  1.3a 

enabled 

3 

92 

Table  6.5  Performance  Improvement  on  CRAY-2  4:1  Abrupt 
Contraction  Flow  (940  elements,  1024  nodes)  Modified  Johnson- 
Segalman  fluids  with  De=0.3. 


CYBER-205  at  U.  of  Minnesota 


EJLUCQDE 

fie 

Comoiler 

Uector 

Iters. 

CPUfsec.) 

O 

O 

0.3 

ftn  200 

disabled 

3 

487 

0.1 

0.3 

VAST  ftn 

enabled 

3 

328 

0.0 

3.0 

VAST  ftn 

enabled 

— 

721 

Table  6.6  Performance  on  CYBER-205  4:1  Abrupt  Contraction 
Flow  (940  elements,  1024  nodes)  Modified  Johnson-Segalman 
fluids. 


7.  Conclusions  and  Future  Directions 

The  following  positive  conclusions  may  be  drawn  from  the  completed 
research: 


•The  pilot  algorithm  has  been  significantly  enhanced  and 
improved  during  the  funding  period. 

•In  a  modern  computing  environment,  particle  tracking 
stress  calculations  can  be  performed  at  reasonable  cost 
on  realistic  grids. 

•The  new  algorithm  has  been  used  to  gain  physical  insight 
in  several  important  flow  problems  relevant  to  polymer 
processing  and  measurement. 
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•Continuing  investigation  of  physical  problems  for  which 
experimental  results  exist  promises  to  shed  further  light 
on  the  modelling  of  polymeric  systems. 

•The  stress  calculator  can  operate  independently  of  the 
pseudotime/Galerkin  method  to  calculate  non-Newtonian 
extra  stresses  from  experimentally  determined  velocity 
fields  or  to  be  embedded  in  other  solution  schemes  for 
the  equations  of  motion. 

On  the  other  hand,  the  techniques  developed  to  date  are  subject  to  the 
following  limitations: 

•The  techniques  are  limited  to  planar  and  (in  the  future) 
axisymmetric  flows. 

•Flows  that  develop  singularities  as  time  evolves  seem 
to  be  inaccessible  to  the  pseudotime/Galerkin  method  for 
single  integral  constitutive  equations. 

•Transient  analysis  using  integral  models  seems 
prohibitively  expensive  in  the  current  computing 
environment. 

However,  a  major  finding  of  the  completed  research  is: 

•The  particle  tracking  technique  can  be  the  core  of  fully 
dynamic  methods  for  differential  models  that  show 
promise  of  being  robust  to  the  development  of 
singularities. 

This  last  point  is  the  focus  of  continuing  AFOSR-sponsored  research.  A 
preliminary  description  of  the  use  of  the  stress  calculator  for  differential 
models  can  be  found  in  Ref.  {10}.  The  fundamental  idea  is  to  extend  the 
techniques  for  dynamic  shear  flows  {5,6,9,11,14}  to  planar  flows  by 
replacing  the  convected  derivatives  with  a  ‘Lagrangian  difference 
quotient’  along  streamlines.  This  combines  Galerkin  finite  element 
techniques  with  a  ‘method  of  characteristics’  solution  for  the  stress. 
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8.  Publications 


Journal  and  Reviewed  Proceeding  Articles  Appearing  or  Accepted  During 

Reporting  Period 

1.  D.  S.  Malkus,  “Finite  element  methods  for  viscoelastic  flow,” 
Viscoelasticity  and  Rheology  (Madison,  1985),  Proceedings  of 
Mathematics  Research  Center  Symposium,  October  1984,  A.  S.  Lodge,  J. 
A.  Nohel  and  M.  Renardy,  co-editors,  Academic  Press,  Inc.,  New  York 
(1985),  391-419. 

2.  with  M.  F.  Webster,  “On  the  accuracy  of  finite  element  and  finite 
difference  predictions  of  non-Newtonian  slot  pressures  for  a  Maxwell 
fluid,”  J.  Non-Newtonian  Fluid  Mechanics.  25,  pp  .93-127  (1987). 

3.  With  M.  E.  Piesha  and  M.-R.  Liu,  "Reversed  stability  conditions  in 
transient  finite  element  analysis,”  Comp.  Meths.  Appl.  Mechs.  Eng. 
68  pp.  97-114  (1987). 

4.  With  X.  Qiu,  “Divisor  structure  of  finite  element  eigenproblems  arising 
form  negative  and  zero  masses,”  Comp.  Meths.  Appl.  Mechs.  Eng. 
66,  pp.  365-368  (1988). 

5.  With  R.  W.  Kolkka,  M.  G.  Hansen,  G.  R.  lerley,  and  R.  A.  Worthing,  “Spurt 
phenomena  of  the  Johnson-Segalman  fluid  and  related  models,”  J.  Non- 
Newtonian  Fluid  Mech.  29  pp.  303-335  (1988). 

6.  With  J.  Nohel  and  B.  Plohr,  “Time-Dependent  shear  flow  of  a  non- 
Newtonian  fluid,”  in  Current  Problems  in  Hyperbolic  Systems:  Riemann 
Problems  and  Computations  (Bowdoin,  1988),  Contemporary 
Mathematics,  B.  Lundquist,  editor,  Amer.  Math.  Soc.,  Providence,  Rl,  to 
appear, 1 989. 

Other  Proceedings  Articles  Appearing  or  Accepted  During  Reporting  Period 

7.  “A  fast  algorithm  for  non-Newtonian  flow,"  in  Transactions  of  the 
Third  Army  Conference  on  Applied  Mathematics  and  Computing,  ARO 
Report  #86-1,  pp.  651-660,  (Atlanta,  1986). 
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8.  A  preliminary  version  of  2.  (above),  in  Proceedings  of  the  Workshop  on 
Computational  Fluid  Dynamics  and  Reacting  Gas  Flows  (Minneapolis, 
1986),  Institute  for  Math,  and  Its  Applies. 

9.  A  short,  preliminary  version  of  6.  (above),  in  Transactions  of  the  Sixth 
Army  Conference  on  Applied  Mathematics  and  Computing,  to  appear 
(Boulder,  1988). 

10.  “New  transient  algorithms  for  non-Newtonian  flows,”  Proceedings  of 
the  7th  Conference  on  finite  Elements  in  Flow  Problems  (Huntsville, 
1989),  to  appear. 


Journal  Articles  Submitted  During  Reporting  Period  (in  Review) 

11.  With  J.  Nohel  and  B.  Plohr,  “Dynamics  of  shear  flow  of  a  non- 
Newtonian  fluid,”  CMS  Technical  Summary  Report  #89-14,  1988 
(submitted  to  J.  Comp.  Phys.). 

Unpublished  Reports  During  Reporting  Period 

12.  With  Minwu  Yao,  “On  hole-pressures  in  plane  Poiseuille  flow  over 
transverse  slots,”  MRC  Technical  Summary  Report  No.  2943, 
Mathematics  Research  Center,  University  of  Wisconsin-Madison 
(1986). 

13.  “Computational  methods  for  viscoelastic  flow,”  uncataloged 
Mathematics  Research  Center  Report,  University  of  Wisconsin- 
Madison  (1987). 

Articles/Reports  in  Progress 

14.  With  J.  Nohel  and  B.  Plohr,  “Phase  plane  and  asymptotic  analysis  of 
spurt  phenomena,”  1989. 

15.  With  Y.-C.  Tsai,  “Stability  analysis  of  implicit-explicit  time 
integration  for  viscoelastic  flow,”  1989. 
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16.  With  M.  W.  Johnson,  “On  the  behavior  of  the  Johnson-Segalman  fluid 
mode!  in  step-strains,"  1989. 

Thesis  by  Sponsored  Research  Assistant 

17.  Following  are  the  the  title  and  abstract  of  the  thesis  research  of  Mr. 
Minwu  Yao,  whose  thesis  research  has  been  supported  by  the  current 
grant.  The  anticipated  completion  date  of  the  thesis  is  Fall,  1989.  The 
abstract  is  in  the  form  approved  by  the  Thesis  Committee  at  the  time 
of  Mr.  Yao’s  preliminary  examination  (May,  1988): 

NUMERICAL  SIMULATION 
OF 

PLANAR  NEWTONIAN  AND  NON-NEWTONIAN  FLOWS 
BY  AN  ENHANCED  PARTICLE-TRACKING  FEM  CODE 

by 

MINWU  YAO 


ABSTRACT 

This  project  concerns  computation  of  numerical  solutions  for 
various  flow  problems  of  viscoelastic  liquids,  using  an 
enhanced  particle-tracking  finite  element  method  (FEM) 
program  -  FLUCODE.  The  capabilities  and  limitations  of  this 
approach  are  illustrated  by  in-depth  analysis  of  several  flow 
problems  of  rheological  interest.  The  code  is  based  on  the 
standard  Galerkin  FEM  formulation  and  includes  numerical 
techniques  that  have  been  developed  for  many  years  by  D.S. 
Malkus  and  his  co-workers.  The  latest  version  of  FLUCODE  is 
now  capable  of  simulating  complicated  non-Newtonian  flows 
and  allows  a  full  range  of  preprocessing  and  postprocessing 
options.  The  purpose  of  the  work  is  of  two-fold:  first  is  to  test 
and  evaluate  the  numerical  techniques  used  in  FLUCODE,  to 
compare  our  numerical  results  with  the  available  experimental 
or  other  computational  data,  so  as  to  benchmark  the 
capabilities  and  limitations  of  the  code  and  the  numerical 
techniques.  The  second  is  to  contribute  to  further 
understanding  of  interesting  physical  problems,  such  as  the 
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2-D  flow  over  a  slot,  journal  bearing  flow  and  contraction 
flow.  The  main  topics  discussed  in  this  preliminary  report 
include: 

•Numerical  observation  on  steady  perturbation  due  to  the 
slot  and  the  excess  pressure  rise  across  the  slot 

•Choice  of  computational  boundary  conditions  for 
perturbed  Couette  flow 

•Non-zero  driving  gradient  on  the  hole  centerline  and  the 
HPBL  prediction 

•Spatial  decay  rate  for  stationary  perturbation  of  Couette 
flow 

•Introduction  to  FLUCODE 

•Incompressible  flows  in  journal  bearings 

•Contraction  flows. 


Books 

18.  With  R.  D.  Cook  and  M.  E.  Plesha,  Concepts  and  Applications  of  Finite 
Element  Analysis,  Third  Edition,  Wiley,  1989. 


9.  Interactions 

Spoken  Papers  Presented  at  Meetings,  Conferences,  Seminars,  etc. 

1.  Talks  at  conferences,  workshops,  and  symposia  corresponding  to  items 
1.,  6.,  7.,  8.,  and  9.  in  the  previous  section.  Item  10.  will  be  presented  in 
April,  1989 

2.  “Recent  progress  in  computations  with  single  integral  models,”  4th 
International  workshop  on  Numerical  Methods  in  Non-Newtonian  Flows, 
Spa,  Belgium,  Jun.,  1985. 

3.  “Zero  and  negative  masses  in  finite  element  vibration  and  transient 
analysis,"  Mathematics  Department  Colloquium,  University  of 
Massachusetts-Amherst,  Amherst,  MA,  Oct.,  1985. 

4.  Two-week  lecture  series  on  "The  finite  element  method  in  structural 
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and  continuum  mechanics,"  Nanjing  Aeronautical  Institute,  Nanjing,  P. 
R.  C.  (led  to  Ref.  {4}),  November,  1985. 

5.  Same  title  as  4.,  Computer  Science/Numerical  Analysis  Seminar, 
University  of  Wisconsin,  Madison,  Wl,  Dec.,  1985. 

6.  “Measurement  relations  for  hole  pressures  and  normal  stress 
differences  in  Rheology,”  Joint  MRC/Applied  Math.  Seminar,  University 
of  Wisconsin,  Madison,  Wl,  Feb.,  1986. 

7.  “Type  change  in  non-Newtonian  fluids,”  Mathematics  Department 
Seminar,  Illinois  Institute  of  Technology,  Nov.,  1986. 

8.  Same  title  as  7.,  Fluids  Research  Oriented  Group  Seminar,  Michigan 
Technological  Institute,  Houghton,  Ml,  May,  1987. 

9.  Numerical  aspects  of  the  ‘High  Weissenberg  Number  Problem,’”  5^ 
International  workshop  on  Numerical  Methods  in  Non-Newtonian  Flows, 
Lake  Arrowhead,  CA,  June,  1987. 

10.  “New  transient  algorithms  for  non-Newtonian  flows,"  Mathematics 
Department  Seminar,  Illinois  Institute  of  Technology,  Nov.,  1987. 

11.  “A  stress  calculator  for  BK2Z  and  similar  single  integral  constitutive 
equations,”  Annual  Meeting  of  the  Society  of  Rheology,  Atlanta,  Feb., 
1988. 


Consultative  and  Advisory  Functions 

1.  The  P.  I.  has  maintained  a  close  advisory  relationship  with  the  non- 
Newtonian  fluids  group  at  Illinois  Institute  of  Technology,  headed  by  B. 
Bernstein.  The  pilot  method  was  developed  jointly  with  Bernstein 
during  the  six  years  the  P.  I.  was  on  the  mathematics  faculty  at  I.  I.  T. 
The  P.  I.  has  been  working  closely  with  E.  T.  Olsen  on  the  development  of 
the  stress  calculator.  The  P.  I.  continues  to  make  regular  visits  to  I.  I. 
T.  for  the  purpose  of  co-ordinating  the  research  effort.  He  served  on  the 
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Ph.  D.  Thesis  committee  of  Dr.  Alan  Altman,  who  recently  received  his 
Ph.  D.  under  E.  T.  Olsen.  His  research  involved  the  development  of  a 
temperature-dependent  material  clock  for  the  stress  calculator. 

2.  The  P.  I.  continues  to  collaborate  closely  with  A.  S.  Lodge  (Dept. 
Engineering  Mechanics,  U.  W.  -  Madison),  who  has  developed  a  laboratory 
measurement  device  for  normal-stress  differences  and  viscosities 
based  on  hole-pressures.  The  collaboration  provides  an  opportunity  for 
the  P.  I.  to  relate  his  numerical  studies  to  experimental  results,  and 
has  provided  Prof.  Lodge  with  independent  verification  of  the  empirical 
measurement  relation  upon  which  his  instrument  relies. 

3.  The  P.  I.  continues  as  a  member  of  the  Executive  Committee  of  the 
Rheology  Research  Center  (U.  W.  -  Madison).  This  provides  a  valuable 
opportunity  to  interact  with  a  broad  spectrum  of  researchers  in  non- 
Newtonian  mechanics  from  Mathematics  Chemistry,  Engineering 
Mechanics,  and  Chemical  Engineering. 
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